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We report on the exact computation of the scaling form of the 1-point function, on the upper-half 
plane, of the height 2 variable in the two-dimensional Abelian sandpile model. By comparing the 
open versus the closed boundary condition, we find that the scaling field associated to the height 
2 is a logarithmic scalar field of scaling dimension 2, belonging to a c — —2 logarithmic conformal 
field theory. This identification is confirmed by numerical simulations and extended to the height 3 
and 4 variables, which exhibit the same scaling form. Using the conformal setting, we make precise 
proposals for the bulk 2-point functions of all height variables. 

I. INTRODUCTION 

Sandpile models are open dynamical systems which generically show a wide variety of critical behaviours 0, • I n 
two dimensions, one may be tempted to analyse the stationary state of a sandpile in terms of a local conformal field 
theory. As most sandpile models have a number of unusual features, the underlying conformal theories are expected 
to be peculiar as well. It is now established that at least some of them fall in the class of logarithmic conformal 
theories (LCFT) 0,0. 

In the model we consider here, which is the model originally introduced by Bak, Tang and Wiesenfeld also 
called the (isotropic) Abelian sandpile model 0, a number of non-trivial features specific to a c = —2 LCFT have 
been successfully compared with exact results 0, HI 03, GH IHL El IO> El EB> El > confirming the relevance of the 
conformal description. 

The Abelian sandpile model is most naturally formulated in terms of height variables. These random variables, 
sitting at the sites of a lattice, take the values 1, 2, 3 and 4, and evolve under a discrete time dynamics. When 
the lattice is finite, there is a unique stationary measure, which is uniform on its support, formed by the so-called 
recurrent configurations 6]. If the continuum limit of this measure is a c = —2 LCFT measure, each of the four 
height variables should converge to a definite conformal field in the scaling limit, in such a way that the multisite 
probabilities are reproduced by conformal correlation functions. 

Early on, fundamental differences among the four height variables emerged. The Bombay trick pj allowed Majumdar 
and Dhar to handle the height 1 variable rather easily, and led to the clear identification of its scaling form with a 
primary field 4>{z,z) of dimensions (1,1) [ll|. Technically the other three height variables are harder. To date, the 
only known result, due to Priezzhev, is an exact integral formula for the 1-site probabilities P2, P3, P4 on the discrete 
plane Z 2 (or infinitely far from boundaries) |18j. As these numbers must be subtracted from the random variables 
before making contact with the conformal fields (which have a zero vacuum expectation value on the plane), they 
give no information on the fields themselves. Even though different arguments have suggested a genuine difference 
between the height 1 and the heights 2, 3 and 4 in the bulk, the exact scaling form of the latters remains an open 
question. 

In this Letter, we provide an answer to this problem. The results which follow show that the height 2 scales like a 
logarithmic conformal field ip(z, %) of dimensions (1,1) — the logarithmic partner of the primary field <fi — , and relying 
on numerical simulations of the lattice model, we argue that this is also the scaling form for the heights 3 and 4. Here 
we content ourselves with giving the main results of the lattice calculations and simulations, and their comparison 
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with the LCFT predictions. The details of the analysis will be presented elsewhere. 
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II. THE HEIGHT 2 ON THE UPPER-HALF PLANE 



Information about the scaling form of the height 2 variable requires the knowledge of some of its correlations. The 
1-site probability Pziz) on the upper-half plane is a disguised correlator of the height 2 with itself (its mirror image), 
and gives the simplest way to assess its scaling behaviour. It can be computed in essentially the same way the 1-site 
probability P2 = linxj-^oc P2(z) has been computed on the full plane by Priezzhev [Tsl ]. 
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The height variables are the most basic microscopic variables, but they are strongly correlated over the whole lattice 
because they must make up recurrent configurations. A useful alternative is to treat the recurrent configurations as a 
global random variable, uniformly distributed, and in turn, the recurrent configurations can be traded for the set of 
spanning trees on the lattice one considers |8(- While the spanning tree picture may be avoided for calculations with 
heights 1, it seems to be unavoidable for the height 2 (and 3 and 4). 

Being interested in height observables, one should translate the height values at certain sites into well-defined 
characteristics of the spanning trees. This can be done generally, but for what follows, we only discuss the case of a 
single height 2. In a spanning tree, any site belonging to the branch(es) grown from a site i is called a predecessor 
of i. It can then be shown 18| that a given recurrent configuration has a height 2 at site i if and only if, in the 
corresponding spanning tree, the site i has exactly one predecessor among its nearest neighbours. The distribution 
on spanning trees being uniform, one has to count the spanning trees which satisfy this (non-local) characterization 
at i. 

Priezzhev has devised a method to compute this number [Isj . It involves counting the trees with a certain local 
property, and also graphs with some other non-local properties. The former can be done by standard methods of 
graph theory, like the introduction of a defect matrix, the calculation then boiling down to a small determinant. The 
latter calculation requires to enumerate the graphs with a loop (a non-local constraint but easily computed using a 
defect matrix) and also what Priezzhev called 0-graphs, which are graphs with an imposed 0-shaped circuit. Both 
non-local contributions are divergent, but the combination which is relevant to the site probabilities is finite. 

We have used this technique to compute the probability to have a height 2 somewhere on the upper-half plane. 
With respect to the plane, this situation brings a few changes. 

The most obvious change is the presence of the boundary, on which boundary conditions must be prescribed. Two 
natural -and the only ones known- boundary conditions may be imposed along the boundaries, open (Dirichlet) 
and closed (Neumann). Their precise definition is not important for what follows and can be found for instance in 
[l6| . More important is the fact that the operator which converts an open boundary condition into a closed one, or 
vice-versa, has been identified as a c = —2 primary field fi(z) of weights — | |12|. It will be a crucial ingredient in our 
proof that the height 2 variable scales like the logarithmic field ip. Indeed the probability that a site have height 2 
will be computed for both an open and a closed boundary condition, but the two are not independent since they can 
be related to each other by the use of the field fi. This is a non-trivial consistency check that involves the 4-point 
function with two /x's and two ^>'s. 

The counting of 0-graphs is also a little bit more complicated on the upper-half plane than on the full plane. 
The rotational invariance is lost and the translational symmetry is broken in one direction. Moreover the way the 
O-circuit can touch the boundary has to be controlled. These circumstances increase the number of different graph 
configurations that have to be computed separately, and complicate the expressions. 

The expressions we have obtained for the probabilities using the 0-graph technique are rather cumbersome, to 
the extent that reproducing them here is unlikely to be helpful. The local and the loop contributions are simple but 
those coming from the 0-graphs are fairly long. They involve six-fold integrals and a summation over the sites of the 
upper-half plane. The horizontal translational invariance enables one to carry out the infinite sum over the horizontal 
coordinate, which in turn allows to do one of the six integrations (on the full plane, the two infinite summations can 
be done, and then two of the six integrations, leaving a four-fold integral |Tsj). The integrand/summand involves the 
determinants of 4-by-4 matrices, and contain a rather complicated dependence in m, the distance to the boundary 
of the reference site. 

The asymptotic analysis of these expressions is long and technical, but yields the asymptotic form of i- > 2 P ( TO ) an d 
^2 1 ( m )> the probability that a site at a distance m from the boundary has a height equal to 2, when the boundary is 
all open or all closed. We have obtained the following results, 

P 2 op (m)-P 2 = -L(a+l + blogm) + ©(m" 3 log fc ™), (1) 

HI V 2 ' 

P 2 cl (m)-P 2 = — ^(a + 61ogm) +0(m- 3 log fc TO). (2) 
m 

The height 2 probability P 2 on the full plane is known exactly as a multiple integral, whose numerical evaluation yields 
P 2 = 0.1739 [l8|]- The constant a is also given in terms of integrals which we could not to evaluate analytically, but a 
numerical integration gives a — 0.0403. As to the other constant, we found the exact value b = hP\ = (1 — 2/tt), 
with Pi = 0.07363 the probability that a site in the plane has height 1 0. 

These results give a clear signal that the scaling field describing the height 2 variable has scale dimension 2 and has 
a logarithmic component. The natural candidate in the c — — 2 LCFT is the field ip, namely the logarithmic partner 
of the primary field 0, with conformal dimensions (1, 1). 

If ip is indeed the scaling field of the height 2, the above dominant terms for the subtracted probabilities should be 
equal to (ip(z, z*)) op and (ip(z, z*)) c \ with z — z* = 2im. The non-chiral 1-point functions on the upper-half plane 
are also equal to the chiral correlators on the full plane, of two chiral ip inserted at z and z* , provided the coefficients 
defining the most general chiral 2-point function are properly chosen [Tsj . 
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In order to be able to relate the open and the closed boundary conditions, we consider the more general expectation 
value of ip on the upper-half plane, for which the boundary condition on the real axis is all open except on the interval 
I = [zi, Z2] where it is closed. It can be expressed as a chiral 4-point function on the plane as 
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where = 1.18894 normalizes the lattice 2-point functions of the fields p, [l2|]. In the sandpile model, it yields 
the (scaling limit of the) probability that the site z has a height equal to 2 in presence of a boundary which is open 
on R \ / and closed on /. The two lattice results Q and J5J are then easily recovered upon taking z\, z 2 — > and 
— zi, Z2 — > +00 respectively. In these limits, one may invoke the operator product expansion of two /x's, which closes 
on different fields in the two cases, since the remaining boundary condition is different [lq. 
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where Id and Ijv are the identity fields on an open resp. closed boundary and con is the logarithmic partner of Hat. 
Thus we obtain the two expectation values on the upper-half plane as 
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Because the primary field p is degenerate at level 2, the chiral correlator (p(zi)p(z 2 )ip(z3)ip(z4)) satisfies a second- 
order differential equation, but the logarithmic nature of ip complicates the matter. It is the logarithmic partner 
of a primary field (p of dimensions (1,1), into which it transforms under a dilatation, and it is not quasi-primary, 
because it transforms into primary fields p(z,z) and p{z,z) of dimensions (0,1) and (1,0) under special conformal 
transformations [2(j • The operator product expansion of ip with the chiral stress-energy tensor reads 



T{z)ip{w, w) — 



p(w,w) ip(w,w) — ^<fr(w,w) dip(w,w) 
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As a consequence, the differential equation satisfied by (p{zi) p(z 2 )ip(z3)ip(z4 : )) 1S inhomogeneous |3J, where the in- 
homogeneity depends on the four other 4-point functions (p,ppip), (p,p,(fyip), (ppipp) an d (pp/ip(f)). These 4-point 
correlators themselves satisfy inhomogeneous differential equations, and depend on yet four other correlators, 
(pppp), (p,p,<pp), (p,pp4>) and (p,p,<fitfi), which involve primary fields only, and hence may be computed independently, 
as solutions of homogeneous differential equations. 

Altogether, there is a sequence of nine correlators to compute, the last one being (pp/t/jip). The five correlators 
which contain a (f> satisfy a first order differential equation, because the <f> is also degenerate at level 2, so each depends 
on one arbitrary coefficient. The other four satisfy second order equations, and depend on two arbitrary coefficients. 
Thus the general 4-point function (p(zi)p(z2)tp(z3)ip(z4 : )} depends on 13 arbitrary coefficients. 

To fix their values, we note, from the operator product expansion in the limit zi, 22 — ► 0, that none of the nine 
correlators mentioned above can have a logarithmic singularity, and that the resulting function (ip(z, z*)) op must be 
a function of z — z* only. These two facts together fix 10 coefficients. One more coefficient is fixed if one imposes that 
the limit —Zi, 22 —> +00 remains regular, i.e. the logarithmic singularity is absent. As is clear from JBJ, this implies 
that (ip(z, z*)) c \ is related, through the image method, to a chiral field ip satisfying (tp(z)tp(z*)) = 0. The same is 
true for the height 1 variable and its field <f> [lq. 

All this leaves just two free parameters, themselves determined by demanding that the limit z\,z 2 — > of (0) 
reproduces the lattice result JIJ. It leads to our final result for the expectation value of -0 on the UHP with a 
boundary closed on [21,22] and otherwise open, 
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2). One readily checks that in the limit 



where x = Zl2Z34 is a cross-ratio of the four insertion points (23 = 2 

Z13Z24 1 y J 

21, 2 2 — > 0, x — > and y/1 — x — > +1, the previous formula reduces to JIJ. 

To make the connection with the probability of a height 2 in presence of a closed boundary, one takes 
+00. Again x goes to 0, but y/l — x now tends to —1 ^(|. Indeed if we set z\ 



-R and 22 = R. 



zi, 22 — > 
we find that 



1 



_ (R-z)(R+z*) 



{R-z*)(R+z) nas um ^ modulus and draws a complete circle around when R varies from to +00. This 
non-trivial monodromy factor brings a global change of sign in (|SJ) with respect to the previous case (21,22 — » 0), 
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and the non-logarithmic b terms drop out inside the curly brackets (the last term does not contribute in either case) . 
The resulting 1-point function coincides exactly with the lattice result (J2J, thereby confirming in a non-trivial way 
the identification of ip as the scaling field for the bulk height 2 variable in the sandpile model. We will add further 
support in the next section. 

We finish this section by giving the expectation value of tp in the upper-half plane, in presence of a boundary closed 
on R_ and open on R + , a formula we will use in the next section. In this case, Z\ — > — oo, 22 — * 0, and a simple 
calculation from (JSJ yields 

(^*)> cl ,o P = -^^{261og 

III. SIMULATIONS ON A STRIP 

We have so far verified the identification of the height 2 variable with the field ip in the two extreme cases, namely 
in the upper-half plane with a homogeneous boundary condition on the real line, either all open or all closed. It may 
be tested more deeply by seeing how the profile of the height 2 probability varies in a region where the boundary 
condition changes. The formula ® corresponds to such a situation, with half the real line open and the other half 
closed, and would describe the way the probability P%(z) is interpolated between its open boundary value and its 
closed boundary value. 

We could however not carry out the lattice analysis for this case as we did for the homogeneous boundary condition, 
so that the exact lattice result is not available. Thus we turned to numerical simulations of the model to measure the 
lattice profile of Pziz), with the extra advantage that the two other probabilities P^{z) and Pi(z) can be obtained 
with no more effort. 

Clearly the simulations can only be done on a finite array, in our case, a rectangle of base M and height N. We have 
imposed an open boundary condition on the bottom boundary, a closed boundary condition on the top boundary, and 
we have evaluated the four probabilities Pi(n), P2(n), P^{n) and Pa(ji) on the vertical line joining the bottom to the 
top boundary, and located at mid-length (we have included P\{n) for completeness and for comparison purposes). 
For these results to be significant (with respect to the CFT approach), both the height N and the ratio r = M/N 
must be large enough so that the scaling regime is approached, and so that the rectangle can be considered as a 
good approximation of the infinite strip. We found that the values N = 50 and r = 4 were sufficient, by verifying 
that larger values did not bring any noticeable changes in the curves. To speed up the calculations, we chose open 
boundary conditions on the left and right boundaries. 

We have determined the 1— site probabilities Pi(n) on a sample of nearly 17- 10 9 recurrent configurations, generated 
from 1680 random recurrent configurations in the following way. Applied to each of these, the dynamics of the model 
produces a new recurrent configuration at each time step. The first 10 5 ones are not sampled, and from then on, every 
50th configuration is sampled, for a total of 10 7 sampled configurations for each initial configuration. At each point n 
of the line joining the two boundaries, the probabilities are the ratios of the number of occurrences by the size of the 
sample. From these numbers, we should subtract in each case the probability P& in the bulk (far from the boundaries) 
before making the comparison with the field-theoretic results. The most distant sites from the boundaries are the 
central sites, but because the strip is finite, the lattice probabilities at the central sites have not yet attained their 
limit value P&. So we have subtracted an effective probability P£ s , computed by a fitting procedure. The results are 
shown as colour dots on the left plots of Figures 1 and 2. 

Despite the fact that the size of the sample (~ 10 10 ) is desperately small compared to the volume of the phase 
space (~ 10 5045 !), the fluctuations in the data are surprisingly small, and follow a normal distribution to a good 
approximation. Based on an analysis of these, we expect errors less than 2% for Pi(n), less than 1% for Pz{n) and 
Ps{n), and less than 0.5% for Pi(n). The precision naturally increases for more frequent events (the bulk values are 
Pi = 0.0736, P 2 = 0.1739, P 3 = 0.3063, P 4 = 0.4461 [3). 

The results pertaining to the height 2 probability P2 (n) may be directly compared with the probability profile on 
the infinite strip obtained from the LCFT picture, on the basis that the height 2 variable goes to the field if) in the 
scaling limit. For that, the corresponding curve on the upper-half plane is simply transformed to the infinite strip, of 
width L, by the conformal transformation w — — logz. 

By integrating the infinitesimal transformation of ip given by J7J, one finds its finite transformation under z — > w, 

r 1 z"(w) z"(w) ] 

tp nevf (w,w) = \z'(w)\ 2 \i} old (w,w) - -log|z'(w)| 2 <j) id(w,w) + „ , 2 , s p id(w,w) + 2 Pold(u>, w) \. (10) 
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Under the above mapping z — exp (ttw / L) of the upper-half plane onto the infinite strip, and taking expectation 
values, one obtains (ip(w, w)} a tn P in terms of (ip(w, w))\jhp, (<f>(w, to))uHP) (p(w,w))uhp and (p(w, u)))tjhp- The 
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FIG. 1: Subtracted probabilities to find a height 1 (in yellow) or a height 2 (in blue) at a site lying along the short medial 
line of a rectangle M x N, and going from an open boundary (left end of the graph) to a closed boundary (right end). The 
dots correspond to data obtained from numerical simulations, with M = 200 and N — 50, while the solid curves represent the 
conformal predictions for the limit M — +oo. On the left are plotted the raw data and the functions (II 21) and 11311 to which 

they can be directly compared. On the right, the same plots, but multiplied by the factor (-|) 2 S ^ s ^n/L) ' revea l the presence 
or the absence of a logarithmic dependence. 



expectation of ip is given in to which that of <fi is related through a dilatation, whereas the expectations of p and 
p can be shown to vanish, 
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Putting all together, we obtained the following formula for the expectation value of ip on an infinite strip of width 
L, and coordinates w — u + iv, 



7r\ 2 cos(irv/L) 
LJ sm 2 {Trv/L) 
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The open boundary is at v = and the closed one at v = L. 

The height 1 variable scales exactly to the primary field <f>, normalized as in ifTTJl ig (that is why we introduced 
the factor — | in J7J). So one quickly obtains 



7r\ 2 cos(nv/L) 
LJ sin 2 (nv/L) 



(13) 



The two functions of v in (|12|l and (|13[1 have been plotted in the left graph of Figure 1, as solid lines. In both cases, 
we have chosen L = N + ^ , for the following reason. In the continuum, v = L is the position of the closed boundary, 
on which the Neumann condition is enforced. On the lattice, this condition is imposed through a mirror symmetry 
about the line n = N + |, so that it is at that point that the normal derivative vanishes. 

One sees that the conformal curves are very close to the profiles obtained from numerical simulations, adding much 
support to the identification of ip as the scaling field for the height 2 variable. Because N = 50 is relatively small, the 
plots with L — N would be slightly shifted near the closed side, resulting in a visible discrepancy at the right end of 
the plots. The gross features of the two curves are essentially the same and do not reveal any qualitative difference 
between the height 1 and the height 2, namely the logarithmic nature of the latter. In order to make the difference to 
appear more clearly, we have plotted, in the right graph of Figure 1, the same curves and the same data but divided 



by the damping prefactor (j-) 



cos(7rn/L) 
A (irn/L) 



, which tends to make the function vanish precisely for the values of v where 



the logarithm becomes important. One then sees that the plot for the height 1 is flat, while the other one shows 
a non-trivial dependence, and agrees well with the simulations (the discrepancy in the central region is due to the 
difference between the true bulk value of P\ and P2 and the effective values used in the subtraction, a difference which 
is amplified by the factor L 2 /(ir 2 cos (ttv/L)); the same factor also amplifies any disparity between the dots and the 
curve). We note that for the height 2, the curve is asymmetric around the middle point, a consequence of the extra 
factor I in the probability near an open boundary, see (QJ and J5J). 

The graphs in Figure 2 concern the height 3 and height 4 probabilities. As before the dots in the left plot represent 
the subtracted probabilities obtained from numerical simulations. Suspecting that the height 3 and 4 variables scale 
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the same way as the height 2, we have fitted these data with the expressions they would follow if that were true, 
namely the expression (|12|) . In other words, we have assumed that the probability profiles for the heights 3 and 4 
across the finite discrete strip of width N = L — i have the form 

+ baog(^sin(^))}, i = 3,4 (14) 

The fit of the three parameters, performed on the 46 most central values of n, yields the following values in each 
case, 

pf = 0.3063, a 3 = -0.01243, b 3 = 0.03810 , (15) 
pf = 0.4461 , a 4 = -0.04636 , 6 4 = -0.05667 . (16) 

By comparison, a fit of the data for the height 2 yields b 2 = 0.01896 for the parameter controlling the logarithmic 
dependence (the exact value is b = Q = 0.01841), and b\ — —0.000388 for the height 1 (the exact value is 0). 

Using the fitted values of the parameters, the curves lfl4^l are plotted on the left graph of Figure 2, as solid lines. 
The rather large values of b 3 and 64, as well as the good agreement between the solid curves and the data, are a very 
strong indication that the height 3 and 4 variables possess the same logarithmic scaling as the height 2 variable. The 
plots in the right part of Figure 2 bear the same relation to the left plots as those in Figure 1, and exhibit a clear 
logarithmic signature. 

If the height 3 and 4 variables scale like the height 2 variable, their associated fields must be ipi = ctiip + Pi(f>, for 
i = 3,4. Using the above formulae for the expectation values and the fitted values in l|15|) and 116|l . we obtain the 
relations bi — bai and a,i = aoti + b[3i, from which we deduce the values a 3 — 2.07, (3 3 — —5.21, — —3.08, /3 4 = 4.22. 

In summary, we propose the following field assignments for the scaling limit of the four height variables 

height 1 : hi (z, z) = 2) , height 3 : h 3 (z, z) = a 3 ijj{z, z) + /?3</>(z, z) , (17) 

height 2 : hi(z, z) = ip(z, z) , height 4 : h 4 (z, z) — a^ip{z, z) + 04<t>(z, z) . (18) 

The numerical values are consistent with the sum being zero. 



p i(») - Pf = (t) 
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IV. HEIGHT CORRELATIONS IN THE BULK 



The results we obtained in the previous sections can finally be used to make definite statements about the height 
correlations on the plane. As far as we know, and apart from those involving heights 1 exclusively, these correlations 
have never been computed nor measured from simulations. 

The correlation of two distant heights hi, hj, located at z\, zi on the plane, is given by the 3-point function of the 
corresponding fields with the field u>, the logarithmic partner of the identity, inserted at infinity, 

Pij{z l2 ) - PiPj = (h i (zi,zi)h j (z 2 , Z2M00)), |*i 2 | »1, i, j = 1,2,3,4. (19) 



For this one needs the 3-correlators ((jxfiw), (<fiipu) and (ipipco), the last one anyway requiring the other two. 
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The most general form of the 3-point function (tpijju)) depends on 10 arbitrary constants. Placing the u> at infinity 
and requiring that the result be translationally invariant reduces it to 

(V>(l)V(2)a;(oc)) = ^fc + 2Slog|z 12 |+Alog 2 |z 12 |}, (20) 

\ z 12\ ^ > 

but also enforces the following correlators 

(mm) = (mm) = (mm) = o , (21) 
(mm^(oo)) = j^u, (mm^(oo)) = r^{B+Aiog\z 12 \}. (22) 

1^12 1 \ z 12\ L ' 

Interestingly, the constant A governs the long distance behaviour of all 2-site height correlations, and is known 
form the calculation of the 2-site probability for heights 1, A = — Pi/2 0. Granting the field assignments made 
above for the heights 3 and 4, we conclude that the leading terms of the 2-site probabilities on the plane are given 
by 

P n (r)-P?~-^L, P li ( r )-P 1 P i ~-^logr, Pij {r)-P iPj ^-^L\oir, i,j = 2,3,4, (23) 

with ct2 = 1, ct3 = 2.07 and = —3.08. All height variables are anticorrelated, except the height 4 which is positively 
correlated to the other three. These tendencies are precisely those computed for the four height variables on an open 
boundary [To| . 
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